
namespace Tuna {

template<typename Tprec, int Dim>
inline bool CDS<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t G_dx = Gamma / dx;
    static prec_t dx_dt = dx / dt;
    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {
	aE (i) = G_dx;
	aW (i) = G_dx;
	aP (i) = aE (i) + aW (i) +  dx_dt;
	sp (i) = phi_0(i) * dx_dt ;
    }
    applyBoundaryConditions1D();
    return 0;
}

template<typename Tprec, int Dim>
inline bool CDS<Tprec, Dim>::calcCoefficients2D () 
{
    prec_t Gdy_dx = Gamma * dy / dx;
    prec_t Gdx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i) 
	for (int j = bj; j <= ej; ++j) {
	    aE (i,j) = Gdy_dx; 
    	    aW (i,j) = Gdy_dx; 
	    aN (i,j) = Gdx_dy; 
	    aS (i,j) = Gdx_dy; 
	    aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j) + dxy_dt;
	    sp (i,j) = phi_0(i,j) * dxy_dt ;    
	}
    applyBoundaryConditions2D();
    return 0;
}

template<typename Tprec, int Dim>
inline bool CDS<Tprec, Dim>::calcCoefficients3D () 
{
    prec_t dyz = dy * dz, dyz_dx = Gamma * dy * dz / dx;
    prec_t dxz = dx * dz, dxz_dy = Gamma * dx * dz / dy;
    prec_t dxy = dx * dy, dxy_dz = Gamma * dx * dy / dz;
    prec_t dxyz_dt = dx * dy * dz / dt;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; 
    sp = 0.0;
	
    for (int i =  bi; i <= ei; ++i) 
	for (int j = bj; j <= ej; ++j) 
	    for (int k = bk; k <= ek; ++k) 
	    {
		aE (i,j,k) = dyz_dx; 
		aW (i,j,k) = dyz_dx; 
		aN (i,j,k) = dxz_dy; 
		aS (i,j,k) = dxz_dy;
		aF (i,j,k) = dxy_dz; 
		aB (i,j,k) = dxy_dz; 
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + aN (i,j,k) + 
		    aS (i,j,k) + aF (i,j,k) + aB (i,j,k) + dxyz_dt;    
		sp (i,j,k) = phi_0(i,j,k) * dxyz_dt ; 
	    }
    applyBoundaryConditions3D();
    return 0;
}

} // Tuna namespace















